subroutine CalcFacesVariables()
use Data_cfoam
use Data_Mesh
use DataTime, only: dt, nstep, ttime, dt2
use DataPerfectGas, only: cv, Gamma, TemperaturePRho
use DataBoundarySource, only: BoundSourcePresent, iBoundSurf, CalcBoundSource
use InletBoundaryConditions, only: InletBCon
implicit none
    include 'mpif.h'

    integer :: i,j,k, l,m, j1, iFace
    integer :: ndeltap,bmeshs,bmeshf
    integer :: iProc,statusp(mpi_status_size) 
    real(8) :: dtminx
	real *8 uf(3), ucf(3)  ! local velocities' basis
	real *8 ucsf(3) ! old value of ucf ! s = staroe!
    real(8) :: alfa, NoSlipCoef
    real *8     cZ(3) ! coord of cell centers
    real(8) :: rr

!    real *8 xu(4,16) ! parameters for inflow boundary conditions
    
    real(8) :: R(3), rMod,c, time, Q,dQdT,p0, a
    real(8) :: U0(3), Pres0, Temp0, Rho0
    integer :: finlet1
    
    real(8) :: r1(3), pt1, pr(3), gradP(3), ut1(3),ur(3,3), divU, Err(3), pt,ut(3), dv, ErrPt
    real(8) :: x,y,z, Umean, Lr, Px, CharN(3)

!    ! distribution of boundary inlet parameters 
!    data (xu(1,j),j=1,16)/ & ! radius jean
!                    0, 0.001714543, 0.003851764, 0.007023664,	0.008882004, 0.010368679,&
!					0.012227018,	0.014457029,	0.016315368, 0.018173711, 0.019660382,&
!					0.021518725, 0.023377064, 0.024425507, 0.024955361, 0.025/
!
!!    data (xu(1,j),j= 1,16)/ & ! radius ps90
!!                    0,0.04800719,0.1078494,0.1966626,0.2486961,0.290323,0.3423565,0.4047968,&
!!                    0.4568303,0.5088639,0.5504907,0.6025243,0.6545578,0.6839142,0.6987501,&
!!                    0.7000002/
!
!    data (xu(2,j),j=1,16)/ & ! velocity
!                    406.7245,420.6993,439.4097,443.9626,448.2266,456.2513,432.3959,333.7299&
!                    ,328.3791,332.5684,339.38,349.9492,358.0045,353.8556,299.6601,0/
!    data (xu(3,j),j=1,16)/ & ! temperature
!                    559.1461,578.8375,602.9079,600.7621,598.8358,595.0297,501.4955,266.8332,&
!                    245.5667,240.5564,237.9458,234.3195,231.6581,233.552,244.145,254.3445/
!
!    data (xu(4,j),j=1,16)/ &
!                    31856.41, 31863.25,31734.49,31424.45,31162.23,30900.34,30449.77,29860.01,&
!                    29099.75,	28129.51,27170.73,25788.25,	25038.45,25548.08,23896.37,22519.15/

     do k=1,nInternalFaces
        i=meshFaceOwner(k)
        j=meshFaceNeighbour(k)
        l=opposingOwnerFaceLabel(k)
        m=opposingNeighbourFaceLabel(k)

        ! SK changed here
!            utmp(6,100)=epsilon
!            ntimemax=20000
!            if(nstep.lt.ntimemax) then
!                propor=(nstep-1.0)/(ntimemax-1.0)
!                utmp(6,100)=1.0d0-propor+propor*epsilon
!            end if


                
        if(NormalChars)then
            CharN(1:3) = meshSf(1:3,k)
        else
            CharN(1:3) = meshCellCenters(1:3,i) - meshCellCenters(1:3,j)
        end if

        call usdefn(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),Us(1:3,m),U(1:3,i),Unew(1:3,i),U(1:3,j),Unew(1:3,j) &
               ,Psnew(k),    Ps(k),    Ps(l),    Ps(m),    p(i),    pnew(i),    p(j),    pnew(j)     &
               ,Rhosnew(k),  Rhos(k),  Rhos(l),  Rhos(m),  rho(i),  rhonew(i),  rho(j),  rhonew(j)     &
        !	           ,Tsnew(k),    Ts(k),    Ts(l),    Ts(m),    T(i),    Tnew(i),    T(j),    Tnew(j)     &
!               ,Db(k),Df(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed)
               ,Db(k),Df(k),CharN,dt,dtminx,epsilonSmoothed)

! 			write(2,'(i3,i4,24e13.5, e20.10, 7e13.5, 8e13.5, 5f8.3,  3e13.5)') nstep, k, Usnew(1:3,k),Us(1:3,k),Us(1:3,l),Us(1:3,m),U(1:3,i),Unew(1:3,i),U(1:3,j),Unew(1:3,j) &
!                   ,Psnew(k),    Ps(k),    Ps(l),    Ps(m),    p(i),    pnew(i),    p(j),    pnew(j)     &
!                   ,Rhosnew(k),  Rhos(k),  Rhos(l),  Rhos(m),  rho(i),  rhonew(i),  rho(j),  rhonew(j)     &
!                   ,Db(k),Df(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed
!
!            write(*,*) '@@@ epsilonSmoothed=', epsilonSmoothed
        
        Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
        ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0    

    end do

    include "netInterface.h"

    do j=1,nBoundaryFields
        if (boundaryMeshName(j)(1:4).eq.'proc') cycle ! goto 10
        bmeshs = 1 + boundaryMeshStart(j) ! Phil, 2012-02-28
        bmeshf = boundaryMeshStart(j) + boundaryMeshSize(j)
        
        if (btype(j).eq.6) then ! wall, full-Slip
            do k = bmeshs,bmeshf
                i=meshFaceOwner(k)
                
                Usnew(1:3,k) = U(1:3,i)
                Psnew(k) = P(i)
                Rhosnew(k) = Rho(i)
                Esnew(k) = Enew(i)
            end do
        end if
        
        if (btype(j).eq.4) then ! wall, full-Slip
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                ! SK changed this to a standard wall bcs..

                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)


                ! SK changed here

!                    utmp(6,100)=epsilon
!                    ntimemax=20000
!                    if(nstep.lt.ntimemax) then
!                        propor=(nstep-1.0)/(ntimemax-1.0)
!                        utmp(6,100)=1.0d0-propor+propor*epsilon
!                    end if

                NoSlipCoef = 1.0d0

                call usdef3(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                       ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                       ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
                !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                       ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,NoSlipCoef)

                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
            end do 
        end if

        if (btype(j).eq.1) then ! Outlet
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)

                ! SK changed here
!                utmp(4,100)=0.0
!!____________________________________________________________________________________________!
                if(ViscPartOnly)then
                    x = meshFaceCenters(1,k)
                    y = meshFaceCenters(2,k) !- 0.5d0
                    z = meshFaceCenters(3,k)
                    
                    Umean = 0.01
                    Lr = 0.005
                    
                    Px = - 12d0 * 1.5e-5 / Lr**2 * Umean
                    
                    utmp(5,j) = Pfon + Px * x
                    utmp(1,j) = -0.5d0 / 1.5e-5 * Px * z*(Lr-z)
                    utmp(4,j) = Tfon
                end if
!!                r(1:3) = meshFaceCenters(1:3,k)
!!                rMod = sqrt(sum(r(1:3)*r(1:3)))
!!                c = sqrt(1.4*PFon/RoFon)
!!                time = ttime + dt2
!!                call CalcQ3D(r,time,Q,dQdt,p0,U0)                
!!
!!                Rho0 = RoFon + (P0-PFon)/c**2
!!                Temp0 = p0 / (Rho0 * cv * (Gamma-1d0))
!!                Pres0 = p0
!!____________________________________________________________________________________________!

                if(TotalPresTempGiven)then
                    call CalcOutletBoundCond_TotalPresTempGiven(Usnew(1:3,k),Psnew(k),Rhosnew(k), &
                                    Unew(1:3,i),pnew(i),Rhonew(i),utmp(1,j))
!                    Usnew(1,k) = utmp(1,j)
!                    Usnew(2:3,k) = 0d0
!                    Psnew(k) = pnew(i)
!                    Rhosnew(k) = Rhonew(i)
                else
                    U0(1) = utmp(1,j) ! UFON
                    U0(2) = utmp(2,j) ! VFON
                    U0(3) = utmp(3,j) ! WFON
                    Temp0 = utmp(4,j)
                    Pres0 = utmp(5,j)

!***********************************************************************************!
!                    r(1:3) = meshFaceCenters(1:3,k)
!!!                    call CalcQ2D_Hankel(r,ttime+dt2,Pres0,U0,Rho0)
!                    call CalcQ3D_ConvectedMonopole(r,ttime+dt2,Pres0,U0,Rho0)
!!
!                    Temp0 = Pres0 / (Rho0 * cv * (Gamma-1d0))
!***********************************************************************************!
                    if(sum(U0*meshSf(1:3,k)).lt.0)then
                        finlet1 = 1d0
                    else
                        finlet1 = 0d0
                    end if
                    call usdef1(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                       ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                       ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
                       ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,Pres0,Temp0,U0(1:3),finlet1) !,0.d0) !utmp(4,100)) !,alfa) 
                end if

                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0  

            end do
        end if

        if (btype(j).eq.2) then ! Inlet
        !	call interp(p1,np1,ttime,uin)
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)

                if(ViscPartOnly)then
                    x = meshFaceCenters(1,k)
                    y = meshFaceCenters(2,k) !- 0.5d0
                    z = meshFaceCenters(3,k)
                    
                    Umean = 0.01
                    Lr = 0.005
                    Px = - 12d0 * 1.5e-5 / Lr**2 * Umean
                    
                    utmp(5,j) = Pfon + Px * x
                    utmp(1,j) = -0.5d0 / 1.5e-5 * Px * z*(Lr-z)
                    utmp(4,j) = Tfon

!                    Pin = 1.0001d5
!                    Pout = 1d5
!
!                    x = meshFaceCenters(1,k)
!                    y = meshFaceCenters(2,k) !- 0.5d0
!                    utmp(5,j) = Pin + (Pout-Pin) * x
!!                    a = 0.5 * (Pout-Pin) !/ nu
!                    utmp(1,j) = 0.5 * (Pout-Pin) *(y**2 - 0.25d0)
!                    !Us(1,i) = 0.5 * (Pout-Pin) * (y-1d0) * y / nu
!                    utmp(2,j) = 0d0
!                    utmp(3,j) = 0d0
!                    utmp(4,j) = Tfon
                end if

                U0(1) = utmp(1,j)
                U0(2) = utmp(2,j)
                U0(3) = utmp(3,j)
                Temp0 = utmp(4,j)
                Pres0 = utmp(5,j)

! SK we don't need this interpolation ?
!                if(PS90)then
                if(InletBCon)then
                    cZ(1:3) = meshFaceCenters(1:3,k)
!                    cZ(1:3) = meshFaceCenters(1:3,k)
                    rr=dsqrt(cZ(2)**2+cZ(3)**2)
!!         if (rr.le.0.7d0) &
!                    if (rr.le.0.025d0) &
!                        call intervt(xu,16,rr,utmp(1,100),utmp(3,100),utmp(2,100))
                        
!                    if (rr.le.0.025d0) then
                        call CalcNozzleBoundCond(cZ,U0,Temp0,Pres0)
!                    end if
                end if

!                utmp(4,100)=1.0d0
!                    utmp(6,100)=epsilon

                if(nstep.lt.nTimeStepAcceleration) then ! ntimemax
                    cZ(1:3) = meshFaceCenters(1:3,k)
                    U0(1) = U0(1)*propor+(1.-propor)*UFON
                    U0(2) = U0(2)*propor+(1.-propor)*VFON
                    U0(3) = U0(3)*propor+(1.-propor)*WFON
                    Pres0 = Pres0*propor+(1.-propor)*PFON
                    Temp0 = Temp0*propor+(1.-propor)*TFON
!                    utmp(1,100)=utmp(1,j)*propor+(1.-propor)*UFON
!                    utmp(2,100)=utmp(5,j)*propor+(1.-propor)*PFON
!                    utmp(3,100)=utmp(4,j)*propor+(1.-propor)*TFON
!!                    else
!!                        utmp(1,100)=utmp(1,100)*cos(0.001*(nstep-ntimemax))
                end if

!____________________________________________________________________________________________!
!                r(1:3) = meshFaceCenters(1:3,k)
!!                rMod = sqrt(sum(r(1:3)*r(1:3)))
!!                c = sqrt(1.4*PFon/RoFon)
!!                time = ttime + dt2
!                call CalcQ2D_Hankel(r,ttime+dt2,Pres0,U0,Rho0)
!!                call CalcQ3D(r,time,Q,dQdt,p0,U0)                
!!
!                Temp0 = Pres0 / (Rho0 * cv * (Gamma-1d0))
!                Pres0 = p0
!                Rho0 = RoFon + (P0-PFon)/c**2
!                Temp0 = p0 / (Rho0 * cv * (Gamma-1d0))
!!!                utmp(1,100) = U0(1)
!____________________________________________________________________________________________!

                if(TotalPresTempGiven)then
                    if(sum(U0*meshSf(1:3,k)).lt.0) then
                        call CalcInletBoundCond_TotalPresTempGiven(Usnew(1:3,k),Unew(1:3,i), &
                                            Psnew(k),Rhosnew(k),meshSf(1:3,k),Pres0,Temp0,U0(1)) ! utmp(2,100),utmp(3,100),utmp(1,100))
                    else
                        call CalcOutletBoundCond_TotalPresTempGiven(Usnew(1:3,k),Psnew(k),Rhosnew(k),Unew(1:3,i),Pnew(i),Rhonew(i),U0(1))
                    end if
                else
                    if(sum(U0*meshSf(1:3,k)).lt.0)then
                        finlet1 = 1d0
                    else
                        finlet1 = 0d0
                    end if
                    call usdef1(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                            ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                            !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                            ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
                            ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,Pres0,Temp0, & ! utmp(1,100),utmp(4,100)) !,alfa)  
                             U0(1:3),finlet1) ! 1d0) ! utmp(4,100)) !,alfa)  
                end if

!                call usdef1_TP(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
!                            ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
!                            !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
!                            ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
!                            ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,utmp(2,100),utmp(3,100),utmp(1,100),utmp(4,100)) !,alfa)  

                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
            end do 

        end if

        if (btype(j).eq.3) then ! Symetry plane ! wall
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)


            ! SK changed here

                NoSlipCoef = 0.0d0

                call usdef3(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                           ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                           ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
            !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                           ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,NoSlipCoef)

                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
            end do 
        end if

        if (btype(j).eq.5) then ! Periodic
            ndeltap=boundaryMeshSize(j)/2
            do k=1+boundaryMeshStart(j),boundaryMeshStart(j)+ndeltap
                i=meshFaceOwner(k)
                j1=meshFaceOwner(k+ndeltap)
                l=opposingOwnerFaceLabel(k)
                m=opposingOwnerFaceLabel(k+ndeltap)
                Df(k)=Db(k+ndeltap)

                ! SK changed here

                call usdefn(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),Us(1:3,m),U(1:3,i),Unew(1:3,i),U(1:3,j1),Unew(1:3,j1) &
                       ,Psnew(k),    Ps(k),    Ps(l),    Ps(m),    p(i),    pnew(i),    p(j1),    pnew(j1)     &
                       ,Rhosnew(k),  Rhos(k),  Rhos(l),  Rhos(m),  rho(i),  rhonew(i),  rho(j1),  rhonew(j1)     &
                !	           ,Tsnew(k),    Ts(k),    Ts(l),    Ts(m),    T(i),    Tnew(i),    T(j1),    Tnew(j1)     &
                       ,Db(k),Df(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed)

                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
!                Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
                Usnew(1:3,k+ndeltap)=Usnew(1:3,k)
                Psnew(k+ndeltap)=Psnew(k)
!                Tsnew(k+ndeltap)=Tsnew(k)
                Rhosnew(k+ndeltap)=Rhosnew(k)
                Esnew(k+ndeltap)=Esnew(k)
            end do
        end if
! 10         continue
    end do
    
    if(BoundSourcePresent)then
        j = iBoundSurf
        bmeshs = 1 + boundaryMeshStart(j) ! Phil, 2012-02-28
        bmeshf = boundaryMeshStart(j) + boundaryMeshSize(j)
        do iFace = bmeshs,bmeshf
            Usnew(3,iFace) = Usnew(3,iFace) + CalcBoundSource(meshFaceCenters(1:3,iFace),ttime+dt)
            a = CalcBoundSource(meshFaceCenters(1:3,iFace),ttime+dt)
            if(abs(a).gt.1d-5)then
                write(*,*) 'GenTurb: ', a,meshFaceCenters(1:3,iFace),ttime+dt
            end if
        end do
    end if
    
    do iFace = 1,nFaces
        Tsnew(iFace) = TemperaturePRho(Psnew(iFace),Rhosnew(iFace))
    end do

end subroutine CalcFacesVariables
